The non-equilibrium quantum state, that is impossible to be captured with a thermodynamic description, is a fascinating topic of current physics research. Rich non-equilibrium physics often stems from the competition of dynamical processes and strong interactions.^[1] It may be encountered in many fields, such as solid-state physics, nonlinear optics, fluids, and even in chemistry and biology.^[2–4] Recent research shows that a quantum system far from equilibrium can exhibit remarkable phenomena that are not possible in equilibrium, e.g., the emergence of multiple steady states and rich exotic phases,^[5–8] the observation of many-body localization,^[9–11] prethermalization,^[12,13] and time crystals.^[14–16]

Recently, studies on Rydberg atoms have attracted a great deal of attention in the field of ultracold quantum gases and condensed matter physics,^[17–21] since Rydberg atoms with a large principal quantum number n have exaggerated atomic properties including long lifetimes and strong interactions,^[22] which have dipole–dipole character at short distances and van der Waals (vdW) character at long distances. The vdW interaction between two atoms in identical Rydberg levels scales as n¹¹ and may lead to a blockade effect for a large enough n:^[23] if one atom has been excited to the Rydberg state, the excitation of the neighboring atoms will be suppressed.^[24–27] The tunable vdW interaction, on one hand, can result in novel phenomena such as supersolid phases,^[28–31] second-order phase transitions,^[32] and nonequilibrium phase transitions;^[8] on the other hand, enables quantum information tasks like the realization of quantum gates,^[33,34] quantum simulators,^[35] quantum entanglement,^[36] sensor,^[37] and single photons.^[38]

In particular, Rydberg atoms trapped in optical lattices can provide a good platform for simulating strongly correlated quantum systems^[21,39] and realizing nonequilibrium quantum phase transitions.^[8] With the development of experimental technologies in ultracold atom manipulation,^[40–44] it is now viable to realize the accurate adjustment of many parameters (like interaction strengths and laser amplitudes or frequencies) on the micro- or nanometer scale, which have motivated more research on the nonequilibrium physics in ultracold atoms.^[45–49] For instance, the efficient trapping of Rydberg atoms in a one-dimensional (1D) optical lattice has been realized in experiments, which extends quantum simulation research to the regime of Rydberg atom physics.^[50] In recent years, the quantum features, e.g., uniform, nonuniform phase, and bistability, of single-species Rydberg atoms in 1D and 2D lattice systems were well investigated both theoretically^[5,51–54] and experimentally.^[55–59] More recently, people further took more realistic factors, such as the effect of the tail of the van der Waals interaction on the nonuniform phase^[60] and the influence of the mixed power-law interactions on bistable^[61] into consider. While the phase diagram of a 1D chain of two-species Rydberg lattices has also been investigated.^[62]

Stimulated by these advances, we consider here a two-dimensional (2D) lattice of two-species Rydberg atoms as shown in Fig. 1, where the trapped atoms along two diagonal lines are excited to the same Rydberg state by different laser fields of Rabi frequencies Ω_A and Ω_B. In the case of Ω_A = Ω_B, our lattice reduces to that in Ref. [5], where the fixed points (steady-state populations) experience stable uniform, unstable uniform, stable nonuniform, or oscillating nonuniform phases along with a bistable region. The uniform, nonuniform, and oscillating phases indicate that the atomic lattice exhibits a uniform Rydberg population everywhere, has two sublattices of different Rydberg populations, and oscillates versus time in Rydberg population, respectively. In the case of Ω_A ≠ = Ω_B, however, the phase diagram becomes quite different with all phases (stable, unstable, or oscillating) being nonuniform so that two islets break away from the mainland of Rydberg excitation. As |Ω_A − Ω_B| increases, the two islets will become smaller and smaller until vanish, while the difference between the fixed points of the two sublattices becomes larger and larger. It is also viable to control the phase diagram by modulating the two laser fields to attain different numbers of the stable and unstable fixed points. Adiabatic dynamics further show that the Rydberg population experiences a bistable behavior only in the mainland if its evolution starts from a stable fixed point in the mainland, while it first evolves in the islets and then jumps to the mainland if its evolution starts from a stable fixed point in the islet.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. Two-species atoms (A-type in blue and B-type in green) are excited to a common Rydberg state from a common ground state by different lasers of Rabi frequencies ΩA and ΩB. When A-type (B-type) atoms are excited to the Rydberg state, B-type (A-type) atoms will experience a shift V in regard of simultaneous Rydberg excitation.

As shown in Fig. 1, we consider a 2D lattice with exactly one two-level atom per site. The atoms along the two diagonal lines are excited to a common Rydberg state from a common ground state by different laser fields of Rabi frequencies Ω_A and Ω_B. Thus, it is convenient to divide the atoms into two species: A-type in blue and B-type in green. The laser–atom coupling scheme is shown on the right side of Fig. 1 where |g〉_A(B) and |r〉_A(B) denote the ground and Rydberg states, respectively, for A(B)-type atoms, and V describes the dipole – dipole interaction between the two Rydberg atoms as a vdW potential. Since the vdW potential decreases rapidly with distance, we consider here only the nearest-neighbor interaction.

The Hamiltonian in the interaction picture under the rotating-wave approximation can be written as (ħ = 1)

where each atom is labeled by (x,y) (x = 0,1,2,…; y = 0,1,2,…), namely, the 2D coordinates of the lattice sites. Δ_(x,y) is the detuning between the laser frequency and the atomic resonance for an atom at position (x,y).

The time evolution of this atomic system is described by the master equation of density operator ρ, 2

where decoherence due to spontaneous emission from the Rydberg state of decay rate γ is described by the Lindblad operator ℒ[ρ]. Under the mean-field approximation, we can neglect the intersite quantum correlations and factorize the density matrix into each site.^[63,64] Then equation (2) reduces to the nonlinear coupling equations

where density-matrix elements ρ_n,rr and ρ_n,gr denote, respectively, the atomic population and coherence at site n = (x,y) with regard to the Rydberg state. The Δ_n,eff = Δ_n = V(∑_{〈n′ ≠ n〉} ρ_n′,rr) is the effective detuning including all contributions arising from the vdW potentials.

Setting ∂ρ_n,rr(gr)/∂t = 0 in Eq. (3), we can obtain the steady-state equations

where lattice sites index n = (x,y) has been replaced with the atom species index, i.e., A (x + y = even) and B (x + y = odd). We have also set Δ_A = Δ_B = Δ for simplicity and the superscript s stands for the steady-state solutions. It is clear that once and are found, and can be easily calculated. By solving the two coupled cubic equations, it is then viable to examine the steady-state solutions as the fixed points of this dimerized atomic system. Note that one should think of the fixed points and as a joint pair because the A-type and B-type atoms are not independent. Moreover, the stability of all fixed points to perturbations of certain parameters should be examined with the Routh–Hurwitz criterion.^[65,66]

In general, the quantum phases determined by the joint fixed points and their stabilities can be classified into three kinds of uniform, nonuniform, and oscillating Rydberg excitation. The uniform phase corresponds to the case that the atoms are uniformly excited over the whole lattice, i.e., . The nonuniform phase corresponds to the case that one sublattice has a higher excitation than the other, i.e., . The oscillating phase corresponds to the case that the Rydberg population always oscillates in time. We will investigate the fixed points of this dimerized system and their stabilities numerically in the next section.

In this section, we numerically investigate the fixed points of the dimerized system, which correspond to the steady-state solutions, and their stabilities. For simplicity, we set γ_A = γ_B = γ and Δ_A = Δ_B = Δ. We consider the weak-interaction case with V comparable to Δ.

In Fig. 2, we show the steady-state solutions and , and their stabilities with the tunable detunings Δ in two different cases: (a) the mono-excitation case of Ω_A = Ω_B, and (b) the dual-excitation case of Ω_A ≠ Ω_B. Note that one should think of the two steady-state solutions and as being a joint pair , since they correspond to the two atom species A-type and B-type.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. Steady-state solutions (blue) and (green) as a function of detuning Δ with V = 5γ and (a) ΩA = ΩB = 1.5γ (green curves coincide with blue curves); (b) ΩA = 1.45γ, ΩB = 1.55γ. Thick solid, thin dashed, and thin solid curves denote stable, unstable, and oscillating fixed points in order. Curves with (without) cross symbols denote the uniform (nonuniform) fixed points.

In the mono-excitation case of Ω_A = Ω_B [see Fig. 2(a)], the system is the same as the one presented by Lee and his coworkers.^[5] We first consider the case (a1) where the detuning is gradually increased from Δ = −3 to Δ = 5. At the beginning, the fixed point starts from a stable uniform phase (thick solid), then undergoes a pitchfork bifurcation at the critical point Δ_c1, in which the stable uniform phase develops into the unstable uniform phase (thin dashed) and the stable nonuniform phase (thick solid). At the critical point Δ_c2, the stable nonuniform phase changes into the oscillatory nonuniform phase (thin solid). At the critical point Δ_c3, the oscillatory nonuniform phase changes back to the stable nonuniform phase (thick solid). At the critical point Δ_c4, there is another pitchfork bifurcation, in which the stable nonuniform phase (thick solid) and the unstable uniform phase (thin dashed) degenerate into the stable uniform phase. At the critical point Δ_c5, the stable uniform phase jumps down to a second stable uniform phase. We then consider the case (a2) where the detuning is gradually decreased from Δ = 5 to Δ = −3. It is clear that the phase transitions in the case (a2) are irreversible as compared to those in the case (a1) because the second stable uniform phase jumps up to the oscillatory nonuniform phase at the critical point Δ_c6. Then we have a bistable region between the critical points Δ_c6 and Δ_c5, which can be further divided into three subregions. In the subregion between points Δ_c5 and Δ_c4, both low and high branches are in the stable uniform phase; in the subregion between points Δ_c4 and Δ_c3, the low branch is in the stable uniform phase while the high branch is in the stable nonuniform phase; in the subregion between points Δ_c3 and Δ_c6, the low branch is in the stable uniform phase while the high branch is in the oscillatory nonuniform phase.

However, in the dual-excitation case of Ω_A ≠ Ω_B, as presented in Fig. 2(b), when a tiny difference between Ω_A and Ω_B is introduced, the phase diagram becomes quite different. All uniform phases (no matter stable, unstable, or oscillatory) become nonuniform and two islets tend to appear. In the mainland, the nonuniform phase (no matter stable, unstable, or oscillatory) , while in the islets, the nonuniform phase (no matter stable, unstable, or oscillatory) . In the bistable region, both high and low branches become nonuniform. Note that, the stable (oscillatory) nonuniform phases in the islets may be not overlapped with those in the mainland. The critical points of the phases in the mainland [marked in Fig. 2(b)] and islets [not marked here] may be different , where the superscript M (I) stands for the mainland (islet), and i = 1,2,3,4. In particular, the stable nonuniform phase in the islets , which implies that the A-type atoms can stably have a higher Rydberg population than the B-type atoms, although the latter suffering from stronger laser driven. It is also worth noting that, the mono-excitation case (a) is an particular case of the dual-excitation case (b), in which the islets and mainland are overlapped when Ω_A = Ω_B.

As shown in Fig. 3, the two islets become smaller and smaller until vanish finally as |Ω_A − Ω_B| increases. In particular, the stable (oscillatory) phases in the islets may be not overlapped with those in the mainland. The difference between the steady populations of A-type and B-type atoms become larger and larger. It is worth noting that, when Ω_B > Ω_A, the islets may keep expanding until overlapped with the mainland as Ω_B −Ω_A decreases to 0. After that, when Ω_B < Ω_A, the steady populations of A- and B-type atoms will suddenly reverse, and the islets will keep shrinking as Ω_A −Ω_B increases (not shown here).

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. Steady-state solutions (blue) and (green) as a function of detuning Δ with V = 5γ and (a) ΩA = 1.25γ, ΩB = 1.75γ; (b) ΩA = 1.2γ, ΩB = 1.8γ; (c) ΩA = 1.1γ, ΩB = 1.9γ; (d) ΩA = 1γ, ΩB = 2γ. Thick solid, thin dashed, and thin solid curves denote stable, unstable, and oscillating fixed points in order.

Now we analyze the number N of steady-state solutions and the number N_S of stable steady-state solutions . Note that one should think of the two steady-state solutions and as being a joint pair , since they correspond to the two atom species A-type and B-type. It is found that (Fig. 3) there may be four regions as the detuning is gradually increased from Δ = −3 to Δ = 5. In the first region, there is one mainland containing one steady-state solution N = 1 (a stable nonuniform phase N_S = 1); in the second region, there are one mainland and two islets containing three steady-state solutions N = 3 (two stable nonuniform phases and an unstable nonuniform phase N_S = 2; or two oscillatory nonuniform phases and an unstable nonuniform phase N_S = 0; or an oscillatory nonuniform phase, a stable nonuniform phase, and an unstable nonuniform phase N_S = 1); in the third region, there is a bistable region containing three steady-state solutions N = 3 (an oscillatory nonuniform phase, an unstable nonuniform phase, and a stable nonuniform phase N_S = 1; or two stable nonuniform phases and an unstable nonuniform phase N_S = 2); in the fourth region, there is one mainland containing one steady-state solution N = 1 (a stable nonuniform phase). As the islets may expand or shrink to vanish, N and N_S may be changed. As the islets shrink, say for example, a region with only one oscillatory nonuniform phase may appear, in which N = 1 and N_S = 0. In particular, as the islets expand, an overlap region between the islets (with the stable nonuniform phase and the unstable nonuniform phase) and the bistable region (with the stable nonuniform high branch and the stable nonuniform low branch and the unstable nonuniform phase) may appear, in which N = 5 and N_S = 3.

Figure 4 shows the number N of steady-state solutions in Δ, Ω_B space, with various fixed Ω_A. As we can see, by tuning the laser Rabi frequencies and detunings, the islet region may expand or shrink, which will give rise to the appearance or disappearance of the overlap region with N = 5.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. A phase diagram on the number N of all steady-state solutions in the (Δ, ΩB) space with (a) ΩA = 0.5γ; (b) ΩA = 1.5γ; (c) ΩA = 2γ; (d) ΩA = 3γ. The dark blue, green, and yellow colors denote N = 1, N = 3, and N = 5, respectively.

Figure 5 shows the number N_S of stable steady-state solutions in (Δ, Ω_B) space, with various fixed Ω_A. As we can see, by tuning the laser Rabi frequencies and detunings, the islet region may expand or shrink, which will give rise to the appearance or disappearance of the overlap region with N_S = 3. These results may have significance for the control of the phase diagram.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. A phase diagram on the number NS of stable steady-state solutions in the (Δ, ΩB) space with (a) ΩA = 0.5γ; (b) ΩA = 1γ; (c) ΩA = 1.5γ; (d) ΩA = 2γ; (e) ΩA = 2.5γ; (f) ΩA = 3γ. The white, blue, light blue, and yellow colors denote NS = 0, NS = 1, NS = 2, and NS = 3, respectively.

In the following, as Δ is adiabatically modulated, we present the dynamical time evolution of the Rydberg steady-state populations in two cases: (1) the evolution starts from a stable fixed point in the mainland [see Figs. 6(a) and 6(c)]; (2) the evolution starts from a stable fixed point in the islet [see Figs. 6(b) and 6(d)]. Note that, to investigate this evolution, the time modulation of the detuning should be adiabatic, so that the Rydberg state populations may evolve to the steady-state ones before next variate δΔ happens. The numerical calculation is based on Eq. (3). The parameters are set as in Fig. 3(b). We first consider the case (1). As the detuning is adiabatically modulated from Δ = −3 to Δ = 5 [solid curves in Fig. 6(c)], the time evolution of starts from a stable fixed point in the mainland [left red points in Fig. 6(a)], goes along the points in the mainland, and then jumps down to the low branch of the bistable points in the mainland at the critical point [solid curves in Fig. 6(a)]. While as the detuning is adiabatically modulated from Δ = 5 to Δ = −3 [dashed curves in Fig. 6(c)], the time evolution of starts from another stable fixed point in the mainland [right red points in Fig. 6(a)], goes along the low branch of the bistable points in the mainland, and then jumps up to the high branch of the bistable points in the mainland at another critical point [dashed curves in Fig. 6(a)]. This is related to the bistable behavior. Note that the steady-state solutions in the black frame in Fig. 6(a) are oscillatory as presented in the inset of Fig. 6(a). Here we use the mean value of the oscillatory solutions as the Rydberg steady-state populations.

Fig. 6.

Figure Option ViewDownloadNew Window

Fig. 6. (a) and (b) Adiabatic evolutions of the steady-state solutions (in blue) and (in green) along with (c) and (d) the common laser detuning Δ attained for the same parameters as those in Fig. 3(b). When Δ is modulated as in panel (c), the adiabatic evolution starts from a stable fixed point in the mainland [see the red points in panel (a)]. When Δ is modulated as in panel (d), the adiabatic evolution of starts from a stable fixed point in the islet [see the red points in panel (b)]. Insets in panels (a) and (b) refer to the oscillating dynamics observed at Δ = 1.33γ and the stable dynamics observed at Δ = 1.33γ, respectively. Solid and dashed curves denote the evolutions versus time T1 and T2, respectively. Black frames denote the regimes where the steady-state solutions are oscillating.

We now investigate the case (2), as shown in Figs. 6(b) and 6(d). With fine-tuned initial population preparations, the time evolution of may start from a stable fixed point in the islet [red points in Fig. 6(b)]. As Δ is adiabatically modulated as the solid curves (or dashed curves) in Fig. 6(d), the time evolution of may go along the fixed points in the islets, and then jump to the fixed points in the mainland. It is worth to note that, when the steady-state solution starts from the fixed points in the mainland, it will not evolve to the fixed points in the islets [compare case (1) with case (2)]. This feature may result in potential applications in the design of novel photonic devices, such as a quantum fuse: when the steady-state population of a system overflows the phase in the islet, the system will exit from the phase and can not reenter it again.

In summary, we have studied the tunable multistability of fixed points in a dimerized system of 2D Rydberg lattice with two laser fields of different Rabi frequencies addressing the trapped atoms along two diagonal lines. It is found that all stable, unstable, and oscillating phases become more and more nonuniform as the difference of the Rabi frequencies gradually increases. This then results in two smaller and smaller islets breaking away from the mainland of steady-state Rydberg populations. Controlled richer phase diagrams with regard to the number of all fixed points and the number of stable fixed points further show that a tristable region may arise in addition to the nonuniform bistable region. We have also examined dynamic evolutions of the Rydberg populations by adiabatically modulating the common laser detuning in two cases corresponding to a starting point in the mainland (1) and in the islets (2), respectively. It is found that the Rydberg population cannot jump out of the mainland in case (1) though exhibiting a bistable behavior, but can jump from the islets to the mainland in case (2) with a reversal jump prohibited. Our findings are expected to be instructive for exploring new nonequilibrium phases in more involved systems of Rydberg lattices.